/*
    Copyright 2007-2011 Patrik Jonsson, sunrise@familjenjonsson.org

    This file is part of Sunrise.

    Sunrise is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 3 of the License, or
    (at your option) any later version.

    Sunrise is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.

*/

/// \file
/// Contains implementation and explicit instantiations of the
/// nonscatter_stage class. \ingroup mcrx

// $Id$

#include "config.h"
#include "mcrx-stage.h"
#include "mcrx-stage-impl.h"
#include "create_grid.h"
#include "arepo_grid.h"
#include "arepo_emission.h"

using namespace std;

bool mcrx::nonscatter_stage::check_stage_state () const
{
  bool state = false;
  try {
    output_file_->pHDU().readKey("MCRX" + stage_ID(), state);
  }
  catch (...) {}
  return state;
}

void mcrx::nonscatter_stage::set_stage_state (bool state)
{
  output_file_->pHDU().addKey("MCRX" + stage_ID(), state, "");
}

void mcrx::nonscatter_stage::setup_objects ()
{
  // First read wavelength vector
  CCfits::ExtHDU& lambda_hdu = open_HDU (*output_file_, "LAMBDA");
  read(lambda_hdu.column("lambda"),lambda_);

  // find reference wavelength
  vector<T_float> templambda(lambda_.begin(), lambda_.end());
  reflambda_= 
    lower_bound (templambda.begin(), templambda.end(),
		 p_.getValue("reference_wavelength", T_float ())) -
    templambda.begin();
  assert(reflambda_>=0);
  assert(reflambda_<lambda_.size());

  cout << "Reference wavelength: " << lambda_ (reflambda_) << " ("
       << reflambda_<< ")\n";

  // Build emission
  cout << "Setting up emission" << endl;
  CCfits::ExtHDU& data_hdu = open_HDU(*input_file_, "PARTICLEDATA");
  emi_.reset(new T_emission(data_hdu,
			    p_.defined("use_reference_for_emission") &&
			    p_.getValue("use_reference_for_emission", bool()))
	     );

  // Build cameras based on -PARAMETERS HDUs in output file
  cout << "Setting up emergence" << endl;
  eme_.reset(new T_emergence (*output_file_));

  // if we are doing kinematic calculation, set the wavelength scale
  // of the cameras
  if(p_.defined("use_kinematics") && 
     p_.getValue("use_kinematics", bool())) {
    T_emergence::iterator stop= eme_->end();
    for (T_emergence::iterator c = eme_->begin();
       c != stop; ++c) {
      (*c)->set_dopplerscale(lambda_);
    }
  }

  // no dust model or grid is used
}


void mcrx::nonscatter_stage::load_file ()
{
  // first set up general stuff
  setup_objects();

  // Only thing to load is camera images
  eme_->load_images(*output_file_, "-NONSCATTER");
}


void mcrx::nonscatter_stage::load_dump (binifstream& dump_file)
{
  // first set up general stuff
  setup_objects();

  // Only thing to load is camera images
  eme_->load_dump(dump_file);
}


void mcrx::nonscatter_stage::save_file ()
{
  // Get normalization and save images
  const T_float normalization = 1./n_rays_;
  eme_->write_images(*output_file_, normalization, m_.units,
		     "-NONSCATTER",
		     p_.defined("compress_images")&&
		     p_.getValue("compress_images", bool ()),
		     p_.defined("write_images_parallel")&&
		     p_.getValue("write_images_parallel", bool ()));

}


void mcrx::nonscatter_stage::save_dump (binofstream& dump_file) const
{
  eme_->write_dump(dump_file);
}


bool mcrx::nonscatter_stage::shoot ()
{
  return shoot_nonscatter ();
}

void mcrx::nonscatter_stage::operator() ()
{
  // this is here so the function is instantiated
  run_stage();
}


